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ABSTRACT 

Newton’s method for nonlinear mechanics problems replaces the governing 
nonlinear equations by an iterative sequence of linear equations. When the 
linear equations are linear differential equations, the equations are usually 
solved by numerical methods. The iterative sequence in Newton’s method can 
exhibit poor convergence properties when the nonlinear problem has multiple 
solutions for a fixed set of parameters, unless the iterative sequences are 
aimed at solving for each solution separately. The theory of the linear 
differential operators is often a better guide for solution strategies in 
applying Newton’s method than the theory of linear algebra associated with 
the numerical analogs of the differential operators. In fact, the theory for 
the differential operators can suggest the choice of numerical linear operators. 
In this paper the method of variation of parameters from the theory of linear 
ordinary differential equations is examined in detail in the context of 
Newton’s method to demonstrate how it might be used as a guide for numerical 
solutions . 


INTRODUCTION 

Nonlinear mechanics problems can be formulated as nonlinear differential 
equations and associated boundary conditions. One approach to solving these 
nonlinear equations is Newton’s method. Newton’s method replaces the nonlinear 
equations with an iterative sequence of linear differential equations. The 
present paper emphasizes that each iteration step consists of two separate 
operations. The first operation, referred to as linearization, is the 
derivation of the linear differential equations. The second operation is the 
solution of the linear equations and is referred to by the name of the method 
of solution for the linear system (e.g., power series, asymptotic series, 
finite-differences, finite-elements, successive approximations, or boundary 
integrals) . 

The emphasis on defining the iteration in Newton’s method as two suc- 
cessive operations is to prevent confusion between Newton’s method and the 
familiar Newton-Raphson method for a set of nonlinear algebraic equations. The 
confusion arises when the second operation is purely numerical and depends on 
a discretization operation. In this case, the operation of discretization 
can be applied to the nonlinear differential operators followed by the linear- 
ization of the Newton-Raphson method. Ortega and Rheinboldt, (ref. 1), prove 
that the operations of linearization and discretization commute. The operations 
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result in the same set of linear algebraic equations for each iteration step 
for both the Newton-Raphson method and Newton f s method. The proof that the 
two operations commute requires that "the discretizations are carried out 
in the same way." 

The present paper examines problems where the discretizations are not 
carried out in the same way. The choice of discrete model is affected by the 
theory of the linear differential equations. Examples of these problems are 
boundary-value problems where multiple solutions exist for a fixed set of 
parameters. This class of problems includes periodic solutions of nonlinear 
dynamics problems and static buckling problems with bifurcation points and 
with limit points. 

Newton’s method, that is, the operation of linearization before dis- 
cretization, supplies two kinds of information for problems with multiple 
solutions. The first kind is qualitative information which is related to the 
convergence of the iterative procedure and is useful in itself. The second 
kind of information is quantitative information that directly affects the 
discrete model. The literature for linear differential equations is vast 
and much of it provides insight into the convergence properties of Newton’s 
method. Rather than attempting a general review of applicable theory, this 
paper examines one method in detail as it relates to Newton’s method. The 
method examined is variation of parameters as it is applied to systems of 
ordinary differential equations. The theory is examined first, followed by a 
discussion of the application of the theory to discrete solutions of nonlinear 
problems with multiple solutions. 

The main body of the paper on variation of parameters is preceded by a 
preliminary section. This section discusses the linearization operation in 
Newton’s method. Once the equations are linearized, different versions of 
Newton’s method receive different names in the literature. These 
versions are briefly reviewed. The section also discussed convergence of 
Newton’s method as it pertains to nonlinear problems with multiple solutions. 
The theoretical results from variation of parameters suggest changes in 
dependent variables that are determined by the given problem and, therefore, 
are applicable to adaptive computer solutions. A final section indicates 
the general nature of such adaptive computer solutions. 


NEWTON’S METHOD 
Fundamental Concepts 

Nonlinear mechanics problems that are formulated as nonlinear ordinary 
or nonlinear partial differential equations can be solved using Newton’s 
method. The basic idea in Newton’s method is to expand the nonlinear operator 
about an assumed or an approximate solution. This expansion yields a new 
nonlinear operator that operates on an unknown correction to the approximate 
solution. It is assumed in Newton's method that nonlinear terms in the 
correction are small compared to linear terms, and the nonlinear terms are 
temporarily neglected. The resulting linear differential equations are 
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solved for an approximate correction which is added to the assumed solution to 
make a new approximation. The procedure is repeated until the corrections are 
small. At each iteration step, the residual error in the solution of the non- 
linear problem is a function of the nonlinear terms neglected in the previous 
iteration step. Convergence of the iterative sequence is almost assured if 
the nonlinear problem has a unique solution. When the nonlinear problem has 
multiple solutions, convergence is not assured in Newton’s method unless 
provisions are made to converge to only one solution for each iteration 
sequence. Examples of problems with multiple solutions are static buckling 
problems with bifurcation points and with limit points and certain nonlinear 
vibrations problems. The theory of linear differential operators is useful 
in guiding numerical computations so that Newton's method converges to the 
desired solution branch. 

The linear operator in Newton’s method is called the Frechet derivative, 
and it is derived from the nonlinear differential operator for the problem. 

Let the nonlinear differential operator be P operating on a scalar function 
or vector function y. The nonlinear problem is 


P(y) = 0 


( 1 ) 


plus associated initial conditions or two-point boundary conditions. Denote 
by y m the approximation to the solution of equation (1) after the mth 
iteration step, and denote by 5ynH-l the correction to y m . Then the Newton 
iteration process solves recursively the equations 


p, » n -i] (S V - - p [Vi ] 


= i + m = 

m m— 1 m 


(2a) 

(2b) 


The operator P’[y m _i] in equation (2a) is the Frechet derivative. Formal 
definitions of the Frechet derivative appear in texts on functional analysis, 
reference (2). For nonlinear differential operators operating on continuous 
functions, the Frechet derivative consists of the linear operators that appear 
in a Taylor series expansion in several variables. The expansion is in terms 
of the dependent variables rather than the independent variables, reference (3). 


Examples of Frechet Derivatives 
For example, consider a single nonlinear equation. 

d 2 y /dy\ 2 

— 2 + Cl — I + X sin y - F sin wt = 0 (3) 

dt \dt / 
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with C, X, F and w constants. The linear variational equation, 
equation (2a) , for equation (3) is 


d 2 6y 


m 


dt 


* “fif) 


d5y 


m 


dt 


+ X (cos y -) <5y = -P[y .] 

m— 1 m m-l 


d 2 y 

P[y 1 ] f 

m_1 dt 2 


1 * frl 


+ X sin y F sin wt 
m-1 


(4a) 


(4b) 


The Frechet derivative for the nonlinear operator is the operator 


d 2 ( ) 


~ + 2C 


dv 


dt 


m-1 

dt 


d 

— ( ) + X cos y , ( ) 
- . m-1 


(5) 


The Taylor series expansion in Newton’s method is readily extended to 
partial differential operators. A second example is a nonlinear strain 
expression 


3u 8w _3u 8u 8v _3v 3w 8w 

£ xz 3z 3x 3z 3z 3x 3z 3x 3z 


( 6 ) 


Then 


[y J(«y ) - 

xz y m-l m 

+ 
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m _j_ m ^ m— 1 

3z 3x 3x 

3u - 36u 3v - 

m-1 m m-1 

3z 3x 3x 


35u 
m 

3z 

36v 3v 36v 

m m-1 m 

3z 3z 3x 


3w - 36w 3w _ 36w 

m-1 n m-1 m 

3x 3z 3z 3x 


(7) 


In addition to the linear variational equations of Newton's method, the 
Frechet derivative appears as part of the chain rule of differentiation. 

If the nonlinear operator is written P(y(x),x) to emphasize that y is a 
function of the independent variable, the total derivative of P is 
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( 8 ) 


The derivative of equation (3) with respect to 


t can be written 


d /dy\ dy d /dy\ dy 

— 2 1 — I + — — ( — ) + ^ (cos y) — = wF cos wt (9) 

dt \dt/ dt dt \dt/ dt 


It is often useful to think of y as a function of parameters in addition to 
the independent variables and the operator as P (y (A,x) ,x, A) . Then 



( 10 ) 


where the dot notation denotes partial differentiation with respect to a 
parameter while holding both the independent and dependent variables fixed. 
If y is considered as a function of t and X is equation (3) , 


d^ /3y\ dy d /3y\ 3y 

— y ( — ) + 2C — — J — I + X cos y — + sin y = 0 (11) 

dt \3A/ dt dt \9A/ dX 


Some versions of Newton’s method make use of equation (10) in solving for 
particular solutions of equations (2a) . 


Different Versions of Newton’s Method 

In this paper, the iterative procedure defined by equations (2) is 
called Newton’s method. Bellman (refs. 4 and 5) gave the procedure the name 
quasilinearization. McGill and Kenneth (ref. 6) use the terminology 
generalized Newton-Raphson operator for the Frechet derivative of nonlinear 
differential operators. The three different names are synonymous for the 
general iterative method. 

When the linear variational equations, equations (2a), are solved by a 
specific algorithm, different writers have coined different names for 
specialized versions of Newton’s method. Perrone and Kao (refs. 7 and 8) 
transform equations (2) to finite difference equations and solve the resulting 
linear algebraic equations by relaxation. This algorithm is called nonlinear 
relaxation by Perrone and Kao. 
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Other versions of Newton T s method are connected with solutions of the 
linear variational equations depending on a parameter. In mechanics, it is 
usual to compute a set of solutions for the nonlinear problems for a given 
set of loads or other parameters. It soon becomes apparent to a user of 
Newton’s method that a solution for one load or parameter is a good zeroth 
approximation for a solution for a nearby parameter (ref. 9). If P(u,A) = 0, 
then u is a good zeroth approximation for the solution for P(y,A + AX) = 0. 
The first iteration of Newton’s method, m = 1 in equations (2), is then 

2 

P'EujCfiyp = -P(u, A + AA) = -[p(u)AA + + ...] (12) 


If the nonlinear operator is linear in the parameter so that quadratic and 
higher order terms in AX do not appear, a particular solution of equation (12) 
follows directly from equation (10). For this case. 


- lx <«> 


y i = u + ax ax (1 4 ) 

If Newton’s method is terminated after one iteration, m = 1, equation (14) 
becomes the zeroth approximation for the next increment on X. Na and 
Turski (ref. 10) call this version of Newton’s method a solution by parameter 
differentiation. When the parameter X is a load parameter and the iteration 
in equations (2) is continued for the iteration counter m > 1, Newton's 
method is also known as the incremental method. Stricklin and Haisler 
(ref. 11) review various versions of this approach when the linear variational 
equations, equations (2a), are solved by the finite element method. 

The idea of continuing known solutions of nonlinear problems into 
nearby neighborhoods gives rise to higher order forms of Newton's method. 

These forms retain higher order derivatives in the Taylor series about a 
known solution. These higher order methods stem from Taylor series expansions 
in the independent variable for Initial value problems. Davis (ref. 12) 
designated this method as a solution by analytic continuation because the 
solutions are capable of extension around singular points in the complex 
domain. Weinitschke (ref. 13) applied a similar approach to solve axisymmetric 
shallow shell equations for particular solutions starting at one boundary. He 
used the Newton-Raphson method to satisfy boundary conditions at a second 
boundary. 

Stricklin, et al. , (refs. 14 and 15) and Noor and Peters (ref. 16) 
combine the idea of analytic continuation with parameter differentiation 
by computing higher partial derivatives with respect to a parameter. 

Stricklin and Haisler (ref. 11) refer to the general scheme as a self-correcting 
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incremental approach to nonlinear problems. Noor does not substitute the 
high partial derivatives in a Taylor series in AA, but uses them as basis 
vectors in a Rayleigh-Ritz solution of the original nonlinear problem. 

The modified Newton's method uses the same Frechet derivative for each 
iteration. Instead of equations (2) , the iterative sequence is 


P , [y C) ]( ( Sy m ) = -p 




(15a) 


y m = 

m 


V-i + 


<Sy. 


m 


m = 1,2,3.., 


(15b) 


Convergence 

The advantage of having the same linear operator for each iteration 
step in the modified Newton's method is offset by slower convergence to the 
solution of the nonlinear problem. Parameter differentiation exhibits 
slow convergence near limit points where 3u/9A is infinite. 

Kantorovich (ref. 2) proves sufficient conditions for the convergence 
of Newton's method. The sufficient conditions are restrictive, but, when they 
apply, the nonlinear problem has a unique solution near the zeroth 
approximation . 

In practical applications of Newton's method, there is not enough 
information to apply Kantorovich ' s convergence criterion. However, it is a 
useful guide because nonlinear problems with unique solutions for a fixed 
range of parameters will exhibit rapid convergence of Newton's method. When 
the solutions are not unique, Newton's method can still converge. The lack 
of uniqueness in the nonlinear solutions is reflected by lack of uniqueness in 
<$y m at some iteration step m. A decision on which solution branch to pursue 
must be made before continuing the iteration. 

The theory of the linear differential operators is a guide for making 
these decisions. Application of the theory of linear ordinary differential 
equations in finding multiple solutions of nonlinear problems is the topic of 
the next section. 


VARIATION OF PARAMETERS 

One method of analysis is examined here in detail to show how linearization 
can influence discretization. The method is variation of parameters which 
is used for finding particular solutions for systems of linear ordinary 
differential equations. The theory for linear ordinary differential equations 
is well understood and is a reliable guide for computing their solutions. The 
variational equations of Newton's method are linear ordinary differential 
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equations when the nonlinear problem is governed by a system of nonlinear 
ordinary differential equations. 

Assume that the nonlinear problem can be written as 

p (y) = ^ + F(y,x,A) =0 a _< x £ b (16) 

plus two-point boundary conditions, 

U(y (a), y (b) , A) = 0 (17) 


The dependent variable y is a vector function with n component 
functions of the independent variable x. There are n boundary conditions 
which may be nonlinear. 

The sequence for the mth iteration step, equations (2), is 


d<5y 


m 


dx 


+ F 'E y m i’ x » A -I 5y m = -p ( y m i» x >^ 

m— 1 m m— l 


(18a) 


U 1 i 

~m-l 


5y (a)' 

m 

i 5y m (b) - 


= -u 


m-1 


(18b) 


y = y + 6y 
m m-1 ra 


ill 1 , 2 , 3 f « I 


(18c) 


The Frechet derivative F' [y m _^,x, A] for this case is the Jacobian 
of the function F(y m _jL,x, A) . The shorthand notation 


Vi ■ p ' [ Vi' x ’ >] 


(19) 


will be used in subsequent equations. 

The solution of the linear variational equations, equation (18a), can 
be written in matrix notation as (see ref. 17, for example) 


Sy = $ C + 
m m~m 


<$ y . 


mp 


( 20 ) 
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The matrix contains n columns, each column contains a linearly 

independent solution (J>(x) of the homogeneous differential equations. 


i£ + j 

dx m-1 


4> = 


0 


( 21 ) 


The vector C m is constant. The vector function fiy— is the particular 
solution of the differential equation system. The method of variation of 
parameters is a method of deriving the particular solution using the solutions 
of the homogeneous equations $ m . The method introduces the change in 
variables 


Sy = $ z 
m mm 


( 22 ) 


where z is a vector function. In the new variables z m , equation (21a) 
is 


dS 
m 

dx 


z + $ 
m m 


dz 
. m 

dx 


+ J n $ z 
m-1 m m 


- P(y m-l’ X ’ X) 


(23) 


The terms multiplying z m vanish identically from equation (23) since each 
column satisfies the homogeneous equations, equation (21). The 

equations 


dz 

* — = -P(y m !,x,X) (24) 

m j m-i 

dx 


are solved by inverting the matrix and integrating each equation. 

Z = c - f S-hc) p(y ,,?,A) dc (25) 

m ~m l m m— i 

*^a 


Substituting equation (25) back into equation (22) completes the 
solution for 6y m . Comparing like terms with equation (20) shows that the 
particular solution is 



= GO 

m 





(26) 
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Whenever the integrals in equation (26) exist, the particular solution 
always exists and depends on both the residual error and the solutions of the 
homogeneous equations, equation (21). 


Boundary Conditions and Compatability 

Determining the values of the constants of integration C m completes 
the mth iteration step. The constants depend on the linearized boundary 
conditions, equation (18b). Substitution of equation (25) into equation (18b) 
yields a set of linear algebraic equations 


B C = D (27) 

~m ~m 


B = U f 
~m ~m 


* (a) 

m 


$ (b) 
m 


(28a) 


D rn = -U T - U ■ \D$y (a); Sy (b)] T (28b) 

~m ~m-l m— 1 mp mp 

Assume that the boundary conditions are posed so that B m is an n x n square 
matrix. Then three possible conditions exist for the solution of equations (27) . 
Let the rank of matrix B m be (n - k) where the number k is the index of 
compatibility (Ince, ref. 18). Let the rank of the matrix obtained by 
augmenting B m with the column vector D m be (n - p) . The three conditions 
are the following: 

1. k = 0 There is a unique solution of equations (27) for the constants 
of integration C m . The mth iteration step is complete with a unique 
correction vector 6y m given by equation (20) . 

2. k < p The algebraic equations are incompatible (also referred to 
as inconsistent) and no solution exists for the constants of integration, 

C . 


3. k = p There are k arbitrary constants in the solutions of 
equations (27) 


B (P = 0 
~m ~m 


j — 1,2, .. k < n 


(29) 
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The solution of the linear variational equations, equation (18), contains k 
arbitrary constants Aj . 


6y = > A. V . + Sy , 

■'m / j j mj ^ml 

j=l 


(30) 


The V . are linear combinations of the solutions of the homogeneous 
equations , 


V . = <f> C 3 
mj m ~m 


j =1,2, . . k < n 


(31) 


The particular solution 5y m ^ is completely determined and satisfies the 
boundary conditions. 


6y ml 


$ C 
mm 


+ 6y. 


mp 


(32) 


The constants of integration C m in equations (32) are solutions of a 
reduced (n - k) square matrix. 

The theory for the compatibility of the boundary conditions is related 
to Kantorovich’s convergence criteria for Newton’s method. The first condition 
k = 0 is part of the sufficient conditions for convergence of Newton’s 
method to a unique solution. Problems with unique solutions converge rapidly 
using Newton’s method and the numerical solution strategy for the zeroth 
approximation and incrementing parameters is not crucial for convergence. 

The second condition k < p indicates a break in the iteration 
sequence. This condition is not usually met in practice. The matrix B m 
becomes ill-conditioned near limit points and approaches the condition k < p, 
but does not satisfy the condition exactly. If the assumed solution 
is modified, the boundary conditions may be shifted to make k = 0 or 
k = p. This modification follows the same lines as the procedure for the 
third condition, k = p, which is considered next. 

When k = p, there are k arbitrary constants of integration Aj in 
the correction function 6y m in equation (30) . If the Aj can be assigned 
values, the iteration can be continued. 
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Multiple Solutions 


If k constants of integration Aj are arbitrary after m iteration 
steps, the nonlinear problem can have multiple solutions for a fixed set of 
input parameters. The procedure for determining the different solutions 
depends on the details of the problem, but the analysis contains general 
guidelines for assigning values for the Aj . The discussion here considers 
the (m+l)th iteration for the case dy^i =0 in equation (30). In other 
words, it is assumed that convergence has been obtained for a solution y m _i 
of the nonlinear problem for a certain value of X = X when the Aj are 
zero. 


The next step is to investigate continuing the iteration with at least 
one of the Aj small but finite. The residual error for the (nri-l)th 
iteration using 5y m from equation (30) is 

P(V*.X> - F"[vi,x.x] (~~j (jf) + -- <33) 


The residual is nonlinear in the constants Aj . The lowest order terms are 
quadratic in the constants unless F"[y m ,x,X] containing second derivatives 
in the Taylor series vanishes. The linear variational equation for the 
(m+l) iteration becomes 


d ^i 

dx 


+ J 

m 


$y. 


m+l 


-P(y m ,x,X) 


(34) 


The Jacobian J m is shifted from J m _i and is also a function of the 
Aj. The difference in the two Jacobians appears if the transformation from 
the mth iteration is applied to the (m+l)th step 


6y = $ z 
y m+l m m+l 


(35) 


Then, equation (34) is transformed to 


dz 


m+l 


dx 




j ,]$ 

m-l J i 


m Z m4-l 



P(y m ,x,A) 


(.36) 
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The difference in the Jacobians is a matrix 

C J - J m •,] = F 'Cy 1 + <$y ,x,X] - F'[y -,,x,X] 
m m-i m— l m m— X 


= F "[y_ + ••• (37) 

m— l m 

The leading terms in an expansion for each coefficient of the matrix is 
linear in the Aj when the residual error expansion starts with quadratic 
terms. The leading exponent is always one less than the leading exponent in 
the residual. 

In theory, equation (36) can be solved by variation of parameters. A 
new set of constants of integration Cm-n appears in the solution to be 
determined by a new set of boundary conditions. 


?m+l Sm+1 


= D 


nH-1 


(38) 


The formal solutions contain the undefined Aj as parameters. In numerical 
solutions, the explicit dependence of the solution on the Aj is not known. 

A procedure that can be implemented for numerical solutions is to partition 
the problem of determining the constants of integration equation (38). 

The partitioning identifies k constants of integration as corrections 
6 Aj on the Aj and partitions equation (38) into a k by k problem to fix 
the Aj and oAj . An iterative procedure is to use trial values for the 
Aj and solve for the 6Aj . The k by k problem is solved when the 6Aj 
vanish for finite real A j . 

The SAj appear as the first k variables of Cm+1 when the v mj 
from equation (30) are arranged in the first k columns of $ m in 
equation (35). Once the Aj and 6Aj are determined, the remaining (n - k) 
equations in equation (38) are a linear set of algebraic equations in the 
remaining (n - k) constants in C m+ -^. 

To avoid converging back to the solution y m -± with the Aj = 0, 
a shift in A is introduced to provide for nonzero, real solutions for 
the constants Aj . 

The shift in A is introduced in the right side of equation (34) 
which is replaced by 

P(y m »x,X) = P(y m ,x,X) + P(y m ,x,X)(X - X) + . . . (39) 
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and the right side of equation (36) is 

-$" 1 P(y m ,x,X) = -<rt>(y ,x,X) - rt(y ,x,X)(A - X) + . . . (40) 

m m ram mm 

Introducing X as a free parameter allows iterating on the reduced k by k 
problem by assigning a value to one Aj and using X as part of the 
minimization on 5ym-KL by making <$Aj zero or small compared to Aj . 

Which Aj to replace as a variable by X depends on the modal coupling 
in the given problem. When k = 1, there is only one Aj - Aj_ so the 
choice is not arbitrary. 

The case k = 1 is a common case for problems with limit points or 
isolated bifurcation points. Therefore, the general iterative procedure 
for k = 1 will be examined further. 


Isolated Bifurcation Point 


For an isolated bifurcation point, y m _i is identified as a solution 
of the nonlinear problem at X = X and k = p = 1. In applications of 
Newton’s method, the solution y m -i can be generated by simply incrementing 
the parameter X in numerical solutions to approach the bifurcation point 
and continue past it. As X is varied in increments in any numerical 
solution using Newton’s method, the boundary condition matrix B approaches 
B^ in equation (29). The rows of B can be rearranged by elementary 
operations so that the coefficients of at least one row, say row k = 1, 
are small as X approaches X. For a bifurcation point, the element k of 
the D vector is identically zero for X near X. 

For a limit point, the kth component of_ D may be small but the small 
divisors of B prevent convergence at X = X without analysis similar to 
the bifurcation analysis described here. 

When y m -i is a solution at the bifurcation point X = X, setting 
Aj_ = 0 so that y m = y m _i and varying X continues the solution so that 


5y 


m+1 


<5y_ _ 

— — (A - A) 
dA 


(41) 


as noted in the discussion on parameter differentiation and equation (13) . 

Having one solution near X, the problem is to investigate the second 
solution of the nonlinear problem whose existence is identified by V-^, 
the nontrivial solution of the linear variational equations. 

A simple procedure for Newton’s method is to assign a value to A^ so 
that y m = y m _i + A^V-l is a known function. Solve the first scalar equation 
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of the set of equations (36) where the equations are arranged so that the first 
element of Zm+1 multiplies the vector function Vp The right side of 
the single differential equation is the first element of equation (40). The 
constant of integration SAi for this equation can be set to zero if 
(A - A) is selected to make the first component of in the boundary 

conditions, equation (38), vanish. 

For numerical solutions, the analysis required by the expansion in 
equation (40) can be circumvented by recognizing that while tbe solution for 
6A^ is nonlinear in Ap it is approximately linear in (A - A) . The 
constant of integration for 5A^ can be_minimized by interpolating linearly 
between solutions for residuals P(y m ,x,A + AA^^) and P(y m ,x,A) where 
AAjq+i is small but arbitrary. This interpolation on A to determine 
A = A^^ for fixed A^ is shown schematically in figure 1. The figure 
represents the surface 6A^ = 6A^(A^,A). It is desired to compute inter- 
sections of the surface and the Aj_ - A plane. The surface is tangent to 
the plane along the A axis, but solutions with A^ = 0 are already assumed 
known by equation (41). Curves on the surface for constant A are nonlinear 
in A^ while curves for constant A^ are assumed to be nearly linear. 
Therefore, one can fix (Ai) m ±0 and approximate A™+l = A m + AA by inter- 
polating on A between the calculated points fiA^L (Aj) m , A m ] and 
6A 1 [(Ai) m , A] to approximate 5A]_[ (A 1 ) m , A^+i] = 0. 

Once Ai is determined and A = \tH -1 th e remaining (n - 1) equations 
of equation (38) can be solved to complete the (m-KL)th iteration. When A^ 
is small, the terms in [j m - J m _^] can be neglected in these equations or 
handled by successive approximations. 

Since [j m - J m _i] is of lower degree than the residual, successive 
approximations should converge rapidly for the remaining (n - 1) equations, 
(ref. 18). 

If the linearized boundary conditions are no longer singular because of 
a finite Ap the usual Newton method iteration can be continued with A 
prescribed as an independent parameter. If the boundary conditions are 
nearly singular, another interpolation on A may be required. 

Once the analysis shows how to start the Newton’s method iteration with 
a finite value of Ap different implicit or explicit numerical integration 
methods can be used to find A as a function of A^. 

NUMERICAL SOLUTIONS 

In applying the results of the theory, numerical solutions do not need to 
follow the analysis in every detail. For example, ’’shooting" methods 
(ref. 19) are used for computing particular solutions rather than variation 
of parameters. When the problem is partitioned into a reduced k by k 
problem, the complete transformation indicated by equation (35) need not be 
carried out. The transformation 


41 



(42) 


Sy 


nri-1 


= w z 


m+1 


is sufficient where w is not singular and contains the k eigenfunctions 
V m j as its first k columns. 

The interpolation on X for fixed can be carried over to nonlinear 

partial differential equations which are solved by matrix methods. A finite- 
element code for general shell problems was used to generate the end-shortening 
u as a function of center deflection w for an isotropic square plate in 
compression. The amplitude w = A^ was held constant at 1.7h where h is 
the thickness while u was varied in the role of X. The interpolation on 
X is indicated by the dashed lines in figure 2. One iteration cycle of 
Newton f s method (m = 1) with input w Q = 1.7h and u = 1.0 u cr (point A 
in the figure) shifted the output to w^ = w Q + 6w^ = 1.32h (point B) . 

Another separate iteration with input at point C shifted the output to 

wj[ = 1.99h and u = 4.47 u C r (point D) . Interpolating on u gives 

Swi = 0 at u = 2.97 u cr (point E) . The latter value of u was sufficiently 

accurate to obtain convergence of Newton's method at point F on the solid 

curve. The solid curve was generated by varying u in increments starting 

at point F using the standard Newton iteration of the computer program. 

Varying the end-shortening with w = 0 over all the plate would merely 
change the in-plane solution which does not couple with the transverse 
equilibrium equation for an initially flat plate. Starting the iteration on 
the computer solution with an assigned amplitude of the lowest buckling mode 
shape is sufficient to achieve convergence on the postbuckled curve. 

When matrix methods are used for the numerical solution of static 
boundary-value problems, a similarity transformation on the tangent stiffness 
matrix is analogous to the change of variables in variations of parameters. 

When the equations are partitioned using only k modes, equations (42), 
the analog is an equivalence transformation in matrix theory. 

The numerical analogies will not be developed in this paper, but the 
general development should parallel the discussion here. Eigenvectors of 
the numerical analog to the Prechet derivative replace eigenfunctions. A 
k by k reduced problem has been suggested by Almroth, et al., (ref. 20). 
Solving the remaining (n - k) equation by successive approximations as 
suggested here preserves linear independence. 


CONCLUDING REMARKS 

Linearization before discretization in Newton's method allows classical 
linear theory to be applied to nonlinear mechanics problems. The linear 
theory provides useful qualitative information that can affect convergence 
of the iterative solution. 

The example in the paper, variation of parameters from the theory of 
linear ordinary equations, shows clearly the interdependence of the residual 
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error and the properties of the linearized operator in Newton* s method. 

Variation of parameters suggests a change of dependent variables for computing • 
particular solutions. The change in variables is also a means of partitioning 
the problem to speed convergence of Newton* s method when the nonlinear 
problem has multiple solutions for a given set of parameters. 

The change of variables can be extended to numerical solutions using 
matrix methods by an equivalence transformation. The new variables arise 
naturally in the problem which is an advantage for writing computer codes 
for matrix solutions. No prior quantitative information on choice of 
variables is required by the program user. 
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